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ABSTRACT 



Context, In the binary system LS I +61° 303 the peak flux density of the radio outburst, which is related to the orbital period of 
26.4960 ± 0.0028 d, exibits a modulation of 1667±8 d. The radio emission at high spatial resolution appears structured in a precessing 
jet with a precessional period of 27-28 d. 

Aims. How close is the precessional period of the radio jet to the orbital period? Any periodicity in the radio emission should be 
revealed by timing analysis. The aim of this work is to establish the accurate value of the precessional period. 

Methods. We analyzed 6.7 years of the Green Bank Interferometer database at 2.2 GHz and 8.3 GHz with the Lomb-Scargle and 
phase dispersion minimization (PDM) methods and performed simulations. 

Results. The periodograms show two periodicities, P { = 26.49±0.07 d (V] = 0.03775 d" 1 ) and P 2 = 26.92±0.07 d (v 2 = 0.03715 d~>). 
Whereas radio outbursts have been known to have nearly orbital occurrence Pi with timing residuals exhibiting a puzzling sawtooth 
pattern, we probe in this paper that they are actually periodical outbursts and that their period is P avcm ge = T^vl ~ 26.70 ± 0.05 d. 
The period P avcra ae as well as the long-term modulation Pbcat = = 1667 ± 393 d result from the beat of the two close periods, the 
orbital Pi and the precessional P 2 periods. 

Conclusions. The precessional period, indicated by the astrometry to be of 27-28 d, is P 2 = 26.92 d. The system LS I +6L303 seems 
to be one more case in astronomy of beat, i.e., a phenomenon occurring when two physical processes create stable variations of nearly 
equal frequencies. The very small difference in frequency creates a long-term variation of period l/(vi - v 2 ). The long-term modulation 
of 1667 d results from the beat of the two close orbital and precessional rates. 

Key words. Radio continuum: stars - X-rays: binaries - X-rays: individual (LS I +61°303 ) - Gamma-rays: stars 



1. Introduction 

The TeV emitting source LS I +61°303 has radio characteris- 
tics that make it unique not only among the small number of 
gamma-ray emitting systems given in the X-ray binary class, 
a class of binary systems where a neutron star or a black hole 
is orbiting around a normal star, but also amo ng the larger 
group of the radio e mitting X-ray binary systems (Fender et al.l 
1 19971 iMassil 120051; iMirabell |2012|) . The peak flux density of 
the radio outburst, which is related to the orbital period of 
26.49 60 + 0.0028 d, exibits a modulation of 1667±8 d dGregorvl 
|2002|) . Some double peaked outbursts when observed at two 
frequencies show different spectral characteristics. There is a 
first outburst with a flat/inverted spectrum and a second opti- 
cally thin outburst associated with different conditions, as in- 
dicated by its high amplitud e, the spectral index, and the Ha 
emission line measurem ents (Mas si & Kaufman Bern add 2009; 
iGrundstrom et alj | 20 07). The complex spectral sequence found 
in LS I +61°303 finds a natural explanation in the disk-jet cou- 
pling model for microquasars: first, there is a continuous outflow 
with a flat or inverted spectrum, th en an event triggers a shock 
in this slow optically thick outflow (Fender et al.ll2004]). and the 
growing shock creates the optic ally thin outburst (Valtaoia et al.l 
119921; lHannikainen et al.l l2006). One of the characteristics that 
make LS I +61°303 unique among the other radio emitting X- 
ray binary systems is that this spectral evolution, between in- 
verted and o ptically thin spectra, may occur twice during the or- 
bital period (Massi & Kaufman Bernado 2009). This agrees well 



with the iBondi & HovTel (fl94i accretion in an eccentric orbit 
(as in LS I +61°303) predicting two events along the orbit as 
shown for LS I +61°303 by several authors (|Tavlor et"aL 1992; 
Marti & Paredesll 19951: iBosch-Ramon et al.ll2006t iRomero etail 



2007). 



The binary system LS I +61°303, for which the nature of the 
compact object has not yet been established (i.e., a black hole 
or a neutron star), shares th e remarkable property of slow radio 
quasi-periodic oscillations dPeracaula et al.l 1 19971 84 min) with 
the tw o black hole microquasars V404 Cvg (|Ha n & Hiellmine 
1 19921 20-120 min) and GRS 1915+105 dPoolev & Fendedl997L 
Rodriguez & Mirabe jfl997l 20^10 min). 

The radio morphology of the system also shows unique char- 
acteristics. The resolved extended structure changes position an- 
gle (i.e., angle of the projection of the jet onto the plane of 
the sky, measured from north through east) with the surpris- 
ing large variation of 60° in only one day dMassi et al.l [2004; 
Dha wan et alj |2006l). Moreover, the jet is sometimes one-sided 
and at other times two-sided. Because of both of these vari- 
ations in position angle and morphology, the hypothesis that 
LS I +61°303 might be a precessing microquasar was brought 
forth (Ma ssi etal.l2 004). The one-sidedness of jets is usually at- 
tributed to relativistic bulk motion along a relatively small angle 
to the line of sight, which leads to Dopp ler boosting of the jet an d 
deboosting of the counterjet emission (lUrrv & Padovanilll995l) . 
A variation of that angle due to precession would cause vari- 
able Doppler deboosting of the counterjet, making it appear for 
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larger angles (double-sided jet) and disappear for smaller values 
of the angle to the line of sight (one- sided jet, blazar like) . 

In 2006, VLBA observations by iDhawan et alj (|2006) mea- 
sured the same la rge rotation of 60°/day in their images as 
Mass i et al .(2004). Some of the VLBA images, s howing again a 
one-si ded structure were, however, interpreted by Dhawa n et alj 
( 2006j) as a cometary tail pointed away from the companion Be 
star, that is in favor of a pulsar model rather than of the pre- 
cessing microquasar model. The reanalysis of this VLBA data 
set and the resulting higher dynamic range of the self-calibrated 
maps ha s actually revealed a double-sided structure in several 
images (Massi et al] |2012|) . Before we illustrate how the re- 
sults in iMassi et aU (120 12l) brought us to the present investiga- 
tion on the precessional period, let us consider two important 
points in the two following paragraphs, the first concerning self- 
calibration and the second the pulsar model. 

Self-ca libration of interferometric data is a well-established 
technique dCornwell & Fomalontll 19991) . It may fail at SNR<4 
(iMartf-Vidal & Marcaidd 120081) or may create spurious sym- 
metrization for unbalanced closure phase triangles (result- 
ing when a very displaced telescope is added to the array) 
(M assi & Aaronl 11999). None of these two cases apply to the 
used set of only VLBA data, where double-sided structures at 8- 
16 cr are present in 6 out the 12 images. Concerning the effects 
on an image of strong va riations of flux density during observa- 
tions ( Ste wart et al.ll201 H) . the source LS I +61°303 shows a low 
radio flux density at all orbital phases, apart from the maximum 
of the long-term modulation. In these epochs a large outburst 
lasting few days occurs around apastron. This means that dur- 
ing the maximum of the long-term modulation one might expect 
a reduction of the dynamic range of the produced maps around 
apastron, i.e., weak feature s will be lost. This is not the case for 
the VLBA observations of IDhawan et al.l §2006) performed to- 
ward the minimum of the long-term modulation. 

We are therefore faced in LS I +61°303 with a changing 
structure from a double-sided structure to a one-sided struc- 
ture. Are suc h variations possibl e also in the pulsar scenario? 
Simulations (Mold 6n"et alJl2012h show that the emission from 
the cometary tail of a pulsar for a particular orientation and incli- 
nation of the orbit, after almost one orbital cycle and at a partic- 
ular orbital phase may look like a double-sided nebula at a fixed 
position angle. In these conditions, the fading and expanding last 
part of the cometary tail may appear detached from the bright- 
est part of the cometary tail which is closer to the orbit. This 
possibility is clearly ruled out for LS I +61°303, where along 
one orbital cycle the VLBA images show a double-sided jet at 
several different o rbital phases and at different position angles 
(IMassi etalj|2012l) . 

We return thus to the scenario of a precessing microquasar, 
where the two-sided jet suffers from variable Doppler boost- 
ing because the precession c ontinuously changes the angle be- 
tween jet and line of sight dMassi et al.l 120041 120121) . Deriving 
the precessional period from the available radio images is not 
straightforward, because the images reflect the variation of the 
projected angle on the sky plane and, therefore, a combination 
of the ejection angle, inclination, and angle of the precession 
cone. The astrometry pro vides less biased results (IDhawan et alj 
2006; IMassi et alJl2012l) . The astrometry of the VLBA obser- 
vations indicates that the peak in consecutive images describes 
a well-defined ellipse, six to s even times larger than the orbit, 
with a period of about 27-28 d (IMassi et al.ll2012l) . Assuming for 
the microquasar in LS I +61°303 the same core-shift effect typ- 
ically observed in blazars, the peak of each image would then 
correspond to that part of the jet where the emission becomes 



optically thick at the obse rving frequency (|Kovalev et al.l l2008). 
Based on this assumption. IMassi et al.l (1201 2l) interpreted the el- 
lipse as the possible cross-section of the precession cone of the 
jet at the distance where the emission at 8.4 GHz becomes op- 
tically thick. The determined time span of 27-28 d to complete 
the ellipse is a first estimate of the precession period. 

The most likely cause for precession of an accretion disk 
of a compact object is an assymetric supernova explosion of 
the progenitor. As a result the compact object could be tilted 
(Fragi le et all l2007h . In this case either the accretion disk is 
coplanar with the compact object and, therefore, subject to the 
gravitational torque of the Be star or, instead, the accretion 
disk is coplanar with the orbit but tilted with respect to the 
compact object which induces, in the context of general rela- 
tivity, Lense-Thirring precess ion if the compact object rotates 
(IMassi & Zimmermannl bOlO). A deep investigation of these or 
other mechanisms of precession requires the knowledge of the 
precession parameters, such as the period of precession and the 
angle of the precession cone. 

In this paper we present a timing analysis of 6.7 years of 
Green Bank Interferometer (GBI) radio data aimed at a more ac- 
curate determination of the precession period. Sections 2. 1 and 
2.2 present the determined period, called Pi. The section illus- 
trates the case that, while the aim of our research was reached, 
i.e., we obtained a more precise value of Pi, we were presented 
with an additional unexpected result: the beating between the or- 
bital period Pi and the precessional period P2 gives rise to a new 
period, P Rvemse = 2/(vi + V2) modulated by l/(vi -V2) = 1667 d, 

1. e., the long-term modulation. Is the found P av erage the period- 
icity of the observed radio outburst? Indeed, whereas in the lit- 
erature Pi is generally referred to as the period of the radio out- 
bursts, it is also well known that there are differences between 
the observed and predicted (for P - Pi) outburst times, and that 
these timing residuals have vs time a sawtooth pattern (Sect. 
2.3). In Sect. 2.4 we present two important results. First, we 
demonstrate mathematically that indeed the sawtooth function 
adjusts Pi to Paveiage- Then we show that the GBI data folded 
with Peerage present an offset of 13 d at the minimum of the 
long-term modulation equal to the sawtooth function and equal 
to that predicted by the beat between Pi and Pi. In the same 
section we also present the corresponding physical scenario. In 
Sect. 3 we discuss the implications of our results for the ob- 
served periodicity in the equivalent width of the Ha emission 
line in LS I +61°303. In Sect. 4 we present our conclusions. 

2. Data analysis and results 

We analyzed here 6.7 years of the NASA/NRAO GBI 
LS I +61°303 database at 2.2 GHz and 8.3 GHz. The 
database covers three periods: 49379.975-50174.710 MJD, 
50410.044-51664.879 MJD, and 51798.333-51823.441 MJD. 
The samples, flux densities at both frequencies with their cor- 
responding errors, in each of the three time intervals are al- 
most continuous with an average of eight observations per 
day. In order to search for possible periodicities we used 
the Lomb-Scargle met hod, which i s very efficien t on irreg- 
ularly sampled data dLombl Il976t IScarglel 119821) . We used 
the algorithms of the UK Starlink software package, PERIOD 
(http://www.starlink.rl.ac.uk/). For both data sets at 
2.2 GHz and 8.3 GHz (Fig. 1) we obtained the same result: 
two periods Pi = 26.49 ± 0.07 d (vi = 0.03775 ± 0.00010 d" 1 ) 
and P 2 = 26.92 + 0.07 d (v 2 = 0.03715 + 0.00010 d" 1 ). The 
statistical significance of a period is calculated in PERIOD 
following the method of Fisher randomization as outlined 
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Fig. 1. Periodograms of 8.3 GHz (filled squares) and 2.2 GHz (circles) data, a: output of Lomb-Scargle method. The Lomb-Scargle 
analysis gives on the y-axis the significance of the frequency. Two frequencies are found at v\ = 0.03775 d~' (P\ = 26.49 d) and at 
V2 = 0.03715 d~' (P2 = 26.92 d). b: output of Lomb-Scargle method for data > 4<x. c: output of PDM. The most likely period yields 
the minimum dispersion and appears as a minimum in the PDM curve, i.e., specular to the maximum in the Lomb-Scargle plot, d: 
output of PDM for data > Act. 



in iLinnel Nemec & Nemed (119851) . The advantage of using a 
Monte-Carlo- or randomization-test is that it is distribution- 
free and that it is not constrained by any specific noise mod- 
els (Poisson, Gaussian, etc.). The fundamental assumption is 
that if there is no periodic signal in the time series data, then 
the measured values are independent of their observation times 
and are likely to have occurred in any other order. One thou- 
sand randomized time series are formed and the periodograms 
calculated. The proportion of permutations that give a peak 
power higher than that of the original time series would then 
provide an estimate of p, the probability that, for a given fre- 
quency window there is no periodic component present in the 
data with this period. A derived period is defined as signifi- 
cant for p < 0.01, and a marginally significant one for 0.01 < 
p < 0.10 (ILinnel Nemec & Nemed 119851) . For both periods 
Pi = 26.49 + 0.07 d (frequency window 0.0374-0.0379 d _1 ) 
andP 2 = 26.92 + 0.07 d (frequency window 0.0369-0.0374d~ I ) 
and for both data sets at 8.3 GHz and 2.2 GHz we obtained 
0.00 < p<0.01. 

2.1. Relative importance between the two periods and 
previous observations 

Figure 1 a shows the results of the Lomb-Scargle analysis for the 
data at 8.3 GHz and 2.2 GHz. There Pi is dominating over P2 
for a factor of 1.8 at 2.2 GHz, and a factor of 1.5 at 8.3 GHz. 
Figure 1 b shows the results of the Lomb-Scargle analysis, if 
only data with flux density > 4<x are used. In this case the two 
periods have a more comparable significance, i.e., there is a fac- 



tor of 1 .4 at 2.2 GHz, and a factor of 1 .2 at 8.3 GHz. We compare 
the Lomb-Scargle results with those ob tained with the phas e dis- 
persion minimization (PDM) method ( Stellingwerf 1978). The 
results of the PDM analysis on the whole data sets are shown 
in Fig. 1 c; the results on data > 4cr are shown in Fig. 1 d. The 
results of the PDM analysis agree very well with those of Lomb- 
Scargle: there is a different significance of the two periods when 
data with low signal-to-noise-ratio (snr) are present in the ana- 
lyzed data set, i.e., Pi dominates. 

This could be the e xplanation why in the pa st the second 
period P 2 was unseen. Taylor & Gregory (1982) found a pe- 
riod of 26.52 + 0.04 d in their radio data set; in 1984 they 
wrote that part of the previous measurements were taken with 
the source in a weak state and repeated the analysis using 
new data and also in cluding the old on es, obtaining the value 
of 26.496 + 0.008 d (ITavlor & Gregorvlll984l) . In 1997 Ray et 
al. reported new observations and gave a period of 26.69 + 
0.02 d, i.e., coincident with our average P aV erage = 7^7; = 
26.70 + 0.05 d, as discussed in Sect. 2.2. These 



0.03775+0.03715 . 

observations (IRav et al.l 1 1997ft are those of our first GBI inter 
val, i.e., 49379.975-50174.710 MJD. In Fig. 2 a one sees that 
these data sample only the inter val of maximum ac t ivity. T he dif- 
ference between the results of iTavlor & Greeorvl (11984 and of 
IRav et all d 1997b is puzzling for both groups. IRav et ai] (1 19971) 
discuss how their best estimate of the period i s significantly dif- 
ferent ( 9cr) from the 26.496 + 0.08 d value of ITavlor & Gre gory 



(ll984 . lGregorv et al.l(11999l), faced with the difference between 
: al.l(ll997l) results, discuss how unlikely 



their value and the Rayet 
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a sudden change in period would be. In the light of our present 
result we see that it is not a sudden change in period but the 
presence of two periods that in the Ray et al. data have compara- 
ble significance. As noted above, the value of iRav etaTI (1997) 
26.69 + 0.02 d corresponds to our P av erage, as discussed in the 
next section. 



2.2. Beating: long-term modulation and P Rvemge 

The two frequencies v x = 0.03775 d" 1 (Pi = 26.49 d) and v 2 = 
0.03715 d" 1 (P 2 = 26.92 d) are only slightly different. This pro- 
duces a beating, i.e., a new frequency is formed v aV erage = ^r 1 , 
modulated with Vb ea t = vi - Vi. For the sum of two sine functions 
the following identity holds 

/ V] — V 2 \ / V] + V 2 \ 

sin (27rvi?)+sin (27rv 2 f) = 2 cos \2n — - — 1\ sin 1 27r — - — 1 1 , (1) 

where the beat frequency (or frequency of the envelope) Vb ea t = 
v\ - v 2 is twice the frequency of the cosine term. In our case, the 

term 7^7 = 0.03775-0.03715 d is eC l Ual t0 1667 ± 393 d Q 

Figure 2 c shows the sum of two sine functions with different 
amplitudes 



/b(f) = sin (2wvi?) + a sin (27TV2?) 

= 2acosl27T — - — 1 1 sin 1 2n — - — 1\ 
+ (1 - a) sin (2nv\f) , 
with a = 0.7, whereas Fig. 2 e shows the function 
fm(t) = (1 + b sin(27rv m ?)) sin(27rvi?), 



(2) 



with b = 0.7, vi = d , and v m 



i 

1667 



(3) 

d _1 . As one can see, 

both Eqs. 2 and 3 are able to reproduce the long-term modula- 
tion; however, the periodograms are rather different. The peri- 
odogram of Eq. 2, shown in Fig. 2 d agrees well with the peri- 
odogram of the GBI data of Fig. 2 b. On the contrary, as one can 
see in Fig. 2 f, in the periodogram of Eq. 3 only Pi = 26.49 d is 
present in the frequency range 0.036-0.039 d -1 . 

2.3. Sawtooth function 



Gregory ( 1999) demonstrated the e xistence of a long-term mod- 
ulation of the peak outburst flux. iGregorv & Neishl (120021) in- 
dicated that the modulation in radio properties may stem from 
periodic ejections of a shell (density enhancement) of gas in the 
equatorial disk of the Be star. 

The long-term modulation is also present in the timing resid- 
uals of the outbursts, i.e., the differe nce between observe d and 
predicted (for P — Pi) outburst time (Greg ory et al.| [l999). The 
observational result is that timing residuals show a surprising 
sawtooth pattern, i.e., with a gradual rise from to about 6 d, 
but then a rapid fall to a large negative value of a bout -7 d. The 
saw to oth function is shown in Figs. 2 and 8 a of Gregor y" 5t al.l 
(1999) and here in Fig. 3 a. In detail the trend is as follows: ob- 
served and predicted outburst times coincide at the peak of the 
long-term modulation (i.e., at the radio maximum) resulting in 
a timing residual equal to r = 0; the timing residual grows lin- 
early with time and at the minimum of the long-term modula- 
tion reaches a maximum of about 6 d where it sharply switches 
to about -7 d. Also surprising is that the transition t =* 6 d to 
t =2 -7d is not related to a strong change in flux density; it 



occurs at a time when the amplitude is low, i.e., at the mini- 
mum of the long-term modulation (Figs. 6 a and 7 b of Gregory 
et al. 1999). After this transition, t starts to grow linearly with 
time reaching the value r = only after about 800 d at the new 
maximum of the long-term modulation. 

We performed a test by using the sawtooth function to cor- 
rect the model of Eq.[3](Pi modulated by 1667 d) and to verify if 
the corrected model is able to reproduce the observed spectrum 
(i.e., Pi and P 2 ). First we have to define the sawt ooth function. 
The sl ope of the sawtooth pattern in Fig. 8 a of IGregorv et al.l 
(119991) . is about 0.008; we generate therefore the sawtooth func- 
tion (Fig. 3 a) with a period of 1667 d as0 



r(t) = 0.008 fmod(r, 1667) 

and include it in Eq.|3] which becomes: 

fm(t) = (l + b sin(27rv m f)) sin(2;rvi(f - r(f))). 



(4) 



(5) 



To study the effect of the sawtooth function without any bias 
resulting from large holes in the sampling, we performed both 
simulations with original sampling and with regular sampling. 
We obtained the same results and here we show those with reg- 
ular sampling. 

The resulting periodogram of the simulated data (Fig. 2 g) is 
shown in Fig. 2h. One sees that the model with the single peri- 
odicity Pi, modulated by 1667 d, once corrected by the sawtooth 
function, is able to reproduce the results of our spectral analysis 
that is the two periods Pi and P%. In the next section we show 
that when one directly uses the two found periods P\ and P%, the 
sawtooth function is naturally explained. 

2.4. Period of the observed outburst and P av erage 

The sawtooth function results from the comparison between ob- 
served and predicted (for P - P\) outburst time. Here we will 
show first analytically (Eq. 6) and then with the GBI observa- 
tions that the observed outburst occurs at P = P aV erage, i-e., at 

1/Vaverage of Eq. 1. 

Analytically, adjusting Pi = 26.49 d (vi = 0.03775 d" 1 ) by 
the given sawtooth function r(f) to 

vi (f - r(f)) = v Y t{\ - 0.008) = (0.03775 x 0.992) t = 0.03745 f(6) 



26.70 d as in Sect. 



0.03745 



gives as a result P aV erage 

2.1, (with variations between 26.65-26.70 d for slopes between 
0.006-0.008). This is an important result. It implies that the tim- 
ing residuals between predicted (at Pi) and observed outbursts 
are equal to the residuals between predicted (at Pi) and P aV erage- 

We therefore ascertained that the observed outburst has peri- 
odicity Paverage- However, if we use Eq. 3 and we simply sub- 
stitute Pi by Paverage, the periodogram fails to reproduce the 
observed periodogram with Pi and P 2 and just shows P aV erage 
(small window on the left in Fig. 2 f). Is P aV erage only an apparent 
periodicity, just produced by Pi and P 2 ? 

Let us fold the GBI data on the periods Pj, P 2 , and P aV erage 
(Figures 4 a-c). The data folded with the orbital period Pi show 
the broad cluster that is well known in the literature (e.g., Fig. 2 c 
in Massi & Kaufman 2009), where flares above 100 mJy occur 
from about phase 0.4 to about phase 0.9. The broadening is due 
to the differences between the observed and predicted (by Pi) 
outburst time that also causes the sawtooth pattern. Figure 4 b, 
shows for the first time the data folded with the precessional 



(7 



) -1 = 1658 + 382 when using > 6 digits 



2 fmod is a function implemented in the math library of C. 
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Fig. 2. Long-term modulation and period analysis, a: 8.3 GHz GBI radio data averaged over 3 d. b: Lomb-Scargle analysis results: 
two frequencies at v\ = 0.03775 d _I (Pi = 26.49 d) and V2 = 0.03715 d _I (P2 = 26.92 d). The small window shows the peak at 
1/1667 d -1 present in all periodograms. c: Sum of two sinusoidal functions at 26.49 d and 26.92 d, with an amplitude ratio 1/0.7. d: 
Lomb-Scargle analysis results: two frequencies at d _I and d l . The significance of the two frequencies becomes identical 
in the periodogram only for an amplitude ratio 1/1. e: Long-term modulation (1667 d) of a 26.49 d periodic outburst, f: Lomb-Scargle 
analysis results: one frequency at 5^49^ ■ The small window to the left shows the peak at P aV erage present in the periodogram of 
a simulation of long-term modulation (1667 d) of a 26.70d periodic outburst, g: Sine wave of periodicity P\, modulated by a sine 
wave of periodicity 1667 d and corrected by a sawtooth function, h: Lomb-Scargle analysis results: two frequencies at d 1 and 
2gij2 ^ 1 as m Figs. 2 b and d. 



period P2. The cluster of the large flares is also evident and su- 
perimposed to scattered smaller flares. If now we fold the data 
with Paverage one would expect, since it is the average of Pi and 
P2, a clustering rather similar to that in Fig. 4 a and b. The re- 
sult, shown in Fig. 4 c, is completely different. First of all, there 
are two clusters. Second, each one of the two clusters is not as 
broad as those with Pi and P2, i.e., the clustering is better. Where 
does the double clustering come from? We used the color green 
for data after 50841 MJD, located in the minimum. A harmonic 
might theoretically give rise to two possible clusters, but in this 
case green and black points had to be present in both clusters. On 
the contrary, all data before the minimum, i.e. the black points, 
cluster at one phase and all points after that, i.e., green points, 
cluster at another phase. The points before and after 50841 MJD 
cluster separately with a shift of 0.5 in phase (or about 13 d). 

We have also folded the simulated data of Eq. 2 with P aV erage- 
In this case the dependency on Pi and P2 is very simple, just two 
sine functions. Nevertheless, the same kind of double clustering 
shown by the GBI data occurs for the simulated data of Eq. 2 
as one can see in the box of Fig. 4 c. The fact that a simple sum 
of sine functions in Pi and P2 produces the same jump, as the 
GBI data, if folded with P aV erage, implies either that this simple 
mathematical form is true in LS I +61°303, or, more likely, that 
the jump is a property of the beating process and its value only 
depends on the two periods Pi and Pi. 

In mathematical terms, as shown in the Appendix, the beat of 
the two sine functions /(Pi) and /(P2) has a phase reset when 
the delay of /(P2) with respect to /(Pi) becomes larger than 
P2/2. Until that point /(Pi) preceeds /(P2) (and P aV erage) giv- 
ing rise to a positive timing residual. After that point /(P2) (and 

1 average 

) preceeds /(Pi). This produces the jump from a large 
timing residual to a nearly equally large, but now negative, tim- 
ing residual observed in the sawtooth function (see Fig. 3 c and 
d and the derived saw tooth function in Fig. 3 b, in impressing 
agreement with the observed one of Fig. 3 a). 



In physical terms the jump is illustrated in Fig. 4 d. The radio 
maximum of the long-term modulation results when the ejec- 
tion, periodical at P = P\, occurs at the smallest angle with 
respect to the line of sight and the Doppler boosting is largest 
dKaufman Bernado et alj I 2002). Because of the precession P2 
the angle of the ejection changes, and the radio minimum corre- 
sponds to an ejection occurring at the largest angle with respect 
to the line of sight. At this point the ejection has travelled half a 
precession cone and turns onto the other half of the precession 
cone, i.e., from I to II in Fig. 4d. This causes the phase jump 
in the data folded with P aV erage and the jump in the saw tooth 
function of Fig. 3 a. 

Finally, in Fig. 4 d, we present the data folded with P aV erage, 
where we fold the data before 50841MJ D (black points) using 
the usual t = 43366.275MJD (lGregorvll2002l) . whereas for the 
data after 50841 MJD (green points) we use to + 13.25 d, correct- 
ing for the jump at the minimum. The better folding of the large 
flares with P aV erage with respect to those with the two real peri- 
odicities Pi and P2 is the observational evidence for the result of 
Eq. 6, i.e., the periodicity of the radio outburst is P aV erage- 

3. Periodicities in the equivalent width of the Ha 
emission line 

Our results have important implications for the short- and long- 
HQ: variations for th e Be star of the LS I +61°303 system. 
IZamanov et ail (1 1999b analyzed Ha spectra of LS I +61°303 and 
determined the same 26.5 d radio period, and in addition deter- 
mined that the Ha emission line equivalent width (EW) varies 
over the s ame time scale as th e long-term radio modulation. 
Moreover. IZamanov et al.l(ll999l) tried to find evidence for long- 
term periodicities in other line parameters like the important B/R 
ratio. 

Most Be/X-ray binaries show asymmetric split Ha profiles, 
the "blue" (B) or "violet" (V) peak and the "red" (R) peak. B/R 
(or V/R) variability refers then to the variation of the relative 
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Fig. 3 . a: Sawtooth func tion based on timing residuals observed 
bv lGregorv et al.l (1 19991) . b: Sawtooth function based on the loss 
of synchronization of the sine functions /(Pi) and f(P<i) of 
Fig. 3 c. c: Sine waves with periods Pi (black line) and P2 (red 
line) starting with synchronized peaks at t = 0. d: Sum of the 
sine waves. The first five bars are at a regular interval equal to 
Paverage, then the next bar follows at a distance of only 13.25 d 
(Pi/2), followed by four bars again at a regular interval equal 
to Paverage as before. The two central bars are symmetric with 
respect to a peak of /(P2) (Fig. 3 c) that is nearly equidistant 
from the preceeding peak and the delayed peak of /(Pi) (see 
Appendix). 



strength of the blue to the red peak. B/R variability cycles are 
common in Be stars forming canonical Be/X-ray binaries con- 
taining accreting X-ray pulsars. As stated above, Zamanov and 
collaborators tried to find evidence in LS I +61°303 for long- 
term periodicities in the B/R ratio, but they conclude that "un- 
fortunately, the resu lts of this search turned out to be negative" 
dZamanovetalJU999l) . 

V/R variability is explained in terms of a nonaxisymmetrical 
equatorial disk in which a one-armed p erturbation (a zone in the 
disk with higher density) propagates (Reig| |20Tl). As a possi- 
ble ex planation for the long-term modulation, Grego ry" & Neishl 



ble expl; 

J2002I) indicated that it may stem from periodi c ejections of a 
shell o f gas in the equatorial disk of the Be star. Gregory & Neish 
(2002) commenting the fact that there are no periodic variations 
in the Ha V/R ratio for LS I +6T303, nevertheless tried to test 
the one-armed density wave model predictions. Their conclu- 
sion is that the radio behavior "is at odds with the predictions of 
a one-armed density wave model. Thus, the one-armed density 
wave model does no t agree quantitatively wi th the measurements 
forLS I +61°303" dGregorv & Neishll2002l) . 

In other words, the origin of the long-term variation of the 
Ha emission line from the disk of the Be star does not seem 
to be related to structural disk variations (no B/R variations), 
whereas it is clearly related to a periodical change in the number 
of emitted photons (the EW(Ha)). A likely candidate as agent of 
these variations is the relativistic precessing jet that could well 
be able to produce EW(Ha) variations of the same timescales as 
those we observe in its radio emission. 




average 



average + offset 



Paveraqe =2/(Vi+ V 2 ) 



#■ # 



Porbital 
1/V! 




Fig. 4. a: Radio light curves (data averaged over 1 d) vs 
^orbit, with a> orbit related to with Pi =26.49 d and 

f =JD2443366.775 dGregorvl 12001 . before (black) and after 
(green) (2400000.5 + 50841 )JD (see Appendix), b: Same as 
related to P2. c: Radio light curves vs 



but for <1> 



precession 



"4 a 

•^average for Paverage = 26.70 d. The small window shows the 
simulated data (Eq. 2) identically folded, d: Radio light curve 
vs <J>average with f changed by Af (At = before (2400000.5 + 
50841)JD and Af = 13.25 d after it, see Appendix), e: Sketch of 
the precessing jet in LS I +61°303 (out of scale). 



4. Conclusions and Discussion 

Our timing analysis and simulations give the following results: 

1. The timing analysis of 6.7 years of GBI data of 
LS I +61°303 results in two frequencies v x = 0.03775 d" 1 
(Pi = 26.49 ± 0.07 d) and v 2 = 0.03715 d" 1 (P 2 = 26.92 + 
0.07 d). The aim of our research, to obtain a better determi- 
nation of the precessional period, indicated by the astrometry 
to be of 27-28 d, has therefore been reached. The period ex- 
ists and it is sligthly above the orbital period. 

2. An additional and totally unexpected result of our research 
is that these two periodicities give rise to a P aV erage = 7-7^ = 
03775+0 03715 = 26.70 ± 0.05 d, modulated with Pb ea t = 



6 



M. Massi and F. Jaron: Orbital and precessional periodicities in LS I +61°303 



— = 0.03775-0.03715 = 1667 ± 393 d ' Ifl ° theI " W01 " dS ' the 

long-term periodicity is equal to the beat of Pi and P2, and 
Pbeat modulates a new period, P aV erage- 

3. We have shown that t he sawtooth function derived in the past 
(Greg ory et alJ [l999) by comparing observed and predicted 
(for P = Pi ) outburst time compares Pi to Paverage- Our result 
therefore confirms the controversial val ue of 26.69 + 0.02d 
obtained in the past by iRav et all (1 19971) for the periodicity 
of the outburst. 

4. Paverage is only an apparent periodicity, a result of the beat 
of Pi and Pi- When the GBI data are folded with P aV erage 
the data, instead of clustering at one specific phase, cluster 
at two phases separated by about 0.5 (13.25 d, i.e., Pi/2) de- 
pending on whether the data are before or after the minimum 
of the long-term modulation. We have shown that the beat of 
Pi and P2 reproduces the same double clustering as the GBI 
data. We find that the time in the minimum of the long-term 
modulation when the phase jump occurs is the time when 
the relative delay between the two functions f(P\ ) and /(P2) 
have reached the maximum delay of P2/2. In a physical sce- 
nario this corresponds, as discussed below, to the point where 
the ejection has travelled half a precession cone and turns on 
the other half of the precession cone. 

5. Paverage and Pb e at are related to each other. They are both pro- 
duced by the beat between Pi and P2, the two real periodic- 
ities. Pbeat, the long-term periodicity, is like P aV eiage, just an 
apparent periodicity. 

6. Pi and the long-term modulation alone, cannot reproduce 
the observed periodogram (i.e., P2). It would be possible to 
reproduce the observed periodogram only by adding a saw- 
tooth function. In this scenario, the sawtooth function must 
be produced by a physical process continuously changing Pi 
to Paverage until a shift in orbital phase of +6 d/26.49 d is 
reached. Then another process would suddenly shift the out- 
burst from +6 d/26.49 d to -7 d/26.49 d in orbital phase, but 
without changing the amplitude of the outburst (remaining 
in the minimum). Instead, we have shown that using directly 
the two periods Pi and P2 that we found in the periodogram 
of the GBI data, one naturally explains the two apparent pe- 
riods Pbeat and Paverage, and the phase jump in P avera ge- 

7. Implications of our results for variations at other wave- 
lengths indicate the precessing relativistic jet as the agent 
responsible for the observed Ho- variations, equal to the vari- 
ations that we observe in radio emission of the jet. 

In conclusion, LS I +61°303 seems to be one more case in 
astronomy of a "beat", i.e., a phenomenon occurring when two 
physical processes create stable variations of nearly equal fre- 
quencies. The very small difference in frequency creates a long- 
term variation of period l/(vi - vi). The first astronomical case 
was that of a class of Cepheids, afterwards called beat Cepheids 
(IOosterhofll957l) . 

The two periods in LS I +61°303 are the precession P2 of an 
ejection having orbital occurrence Pi. Following the results of 
the ra dio spectral index analysis bv lMassi & Kaufman Bernadol 
(2009), the steady ejection of relativistic electrons associated to 
the compact object in LS I +61°303 increases until it termi- 
nates in a transient jet twice along the orbit. As quoted in the 
introduction, this agrees well with the iBondi & HovTel dl944h 
accretion in an eccentric orbit that predicts two events in the 
system LS I +61°303: on e around periastron, and the second 
shifted towards apastron ([Ta ylor et al .l Il992t jMarti & Paredesl 
ll995UBosch-Ramon et al1l2006l:lRomero et alj|2007h . However, 
the ejected relativistic electrons around periastron suffer strong 



inverse Compton losses because they are exposed to stellar pho- 
tons, and only at the second accretion/ejection peak, occuring 
rather displaced from periastron, the relativistic electrons survive 
inverse Compton losses and produ ce the large observed radio 
outburst with period Pi = 26.49 d dBosch-Ramon et al] 2006). 



This ejection periodically changes direction (Ma ssi et al 



2012) 



with a period established in the present paper of 26.92 d. The ra- 
dio maximum of the long-term modulation results when the ejec- 
tion occurs at the smallest angle with re spect to the line of sight 
and t he Doppler boosting is largest dKaufman Bernado et al.l 
2002). The radio minimum results when the ejection occurs at 
large angles with respect to the line of sight. The point where 
the ejection has travelled half a precession cone and turns on the 
other half of the precession cone gives rise to the jump in phase 
for the apparent P aV erage- 

This research has brought up many questions in need of fur- 
ther investigation. Future work needs to be done to establish the 
radio vs Ha relationship, the possible physical processes un- 
der the observed precession, and the Doppler boosting effects. 
Future observations are needed to estimate the angle of the pre- 
cession cone; an estimate of this angle could result from compar- 
ing VLBA observations at the same orbital phase but interlapsed 
(1667/2)d, that is by about 31 orbital cycles. From the different 
position angle of the radio structures one could infer the aperture 
of the cone. Of course the two images interlapsed by 3 1 cycles 
should be done with high precision at the same orbital phase. 
Compare in Fig. 1 in Ma ssi et al. 1 (12011 image A (or image B) 
with image I (or J) taken only one cycle later than A (B), but 
at a sligtly differente phase AO = 0.016. The sma ll phase dif- 
ferenc e already causes differences in the images. In Alber t et al.l 
(120081) it is discussed that two images taken 10 cycles apart and 
at the same orbital phase O = 0.62 are higly similar. Of course 
10 cycles are still far from the requested 31 cycles. However, 
the slightly different position angle between the two images in 
Fig. 1 of those authors indicates that this comparison could be a 
good investigative tool to estimate the precession angle in future 
observations. 
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Appendix A: Beat and sawtooth function 

Let us assume two sine functions f\ and f% with period Pi = 26.49 d and Pi = 
26.92 d. Let us start during the radio maximum, i.e., with zero timing residual 
between the peaks of f\ and fa . In Fig. 4 e a zero timing residual corresponds to 
an ejection with the smallest angle with respect to the line of sight, that is with 
the strongest Doppler boosting. Since Pi is shorter than P2, each cycle /(P2) 
has a delay of P2 - P\ = 0.43 d. The period formed by the beat P av erage has 
a timing residual r between its peak and that of /(Pi) half of that delay, i.e., 
0.43/2 = 0.21 d. This agrees well with r = 0.008 X Pi = 0.21 d of Eq. 4. 

/(/'average ) is always between f\ and fi and its peaks have a regular distance 

of 

^average — 26.70 d (Fig. 3 d). However, at the minimum the accumulated de- 
lay of /(P2) with respect to /(Pi) is almost P2/2. That is, as shown in Fig. 3 c, 
the peak of /(P2) is nearly equidistant from the two peaks of /(Pi), the peak 
about 13 d before and the peak about 13 d after it. The beat with the first peak 
gives rise to a peak of /(Paverage) at about 827.8 d, whereas the beat with the 
second peak gives rise to a peak of /(Paverage) at 841.05 d. The difference be- 
tween these two consecutive peaks of /(Paverage) is 13.25 d. This corresponds to 
about 0.5 in phase, when the data are plotted in phase for Paverage (Fig- 4 c), and 
also corresponds to the jump of about 13 d in the sawtooth function observed by 
iGregorv et al.lfl999T) . After that point, /(P2) is preceeding /(Pi). In the physical 
scenario of Fig. 4 e this jump, or reset of the phase of Paverage, corresponds to the 
point where the ejection has travelled half a precession cone and turns onto the 
other half of the precession cone. Figure 3 b shows r vs time, resulting from this 
simple analysis based on the sine functions of Fig. 3; the resulting slope of the 
function is indeed 0.008, as in Fig. 3 a. 
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